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Introduction 

We  would  like  to  find  a numerical  solution  to  a single  non-linear  1st  order  ODE  (ordinary 
differential  equation)  of  the  form: 
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dy 

dt 


,y,t 


= o 


or: 


dy_ 

dt 


f{y.t)- 


The  basic  philosophy  is  to  start  with  a known  value  of  j/(t)  and  then  approximate  the  next 
value  j/(t  + h ) using  a finite  difference  approximation  to  the  derivative  dy  / dt . For  a 1st 
order  finite  difference: 


dy  _ /n+1)-/n) 

dt  h 


f(yt) 


/n+1)=yM+h-f(y\t*) 


where  yw  refers  to  the  value  of  y from  the  numerical  scheme  at  the  n-th  time  step.  The 
trick  is  to  find  some  appropriate  values  for  h and  f[y  *,t*)  that  lead  to  high  accuracy  & 
stability  without  requiring  an  unmanageable  number  of  function  evaluations  and  time 
steps. 

The  techniques  that  we  will  look  at  will  also  be  applicable  to  systems  of  non-linear  1st  order 
ODEs: 


^=fi{yi>y2>--->yN) for  /=i,2,...,jv 

or  in  vector  notation: 

rfy  = f 

dt 

where  y refers  to  the  set  of  variables  [y1,y2>--->yNf  and  f refers  to  the  set  of  defining 
functions  [/1,/2,...,/JV]  . Note  that  a dependency  on  time  can  always  included  by  defining  a 
new  variable: 

dyN+i  _ i 

dt 

Also,  initial  value  problems  of  order  higher  than  1st  order  can  always  be  recast  as  a system 
of  1st  order  equations  by  defining  new  variables  to  account  for  the  higher  order  derivatives. 
For  example,  the  single  3rd  order  ODE: 
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^-J^-  + 4^-J^-  + 5—  + 2y  = 2sin(t)  with  j/(0)  = j/'(0)  = y"(0)  = 0. 
eft  eft 

Can  be  recast  with  the  following  definitions: 


dy4 

dt 


= 1 


so  the  original  ODE  becomes: 

^ = 2sin(y4)-4y3-5y2-2y1 

with  yq  (0)  = y2  (0)  = y3  (0)  = y4  (0)  = 0 . Now  the  single  ODE  has  become  the  following 
system  of  equations: 


Ti 

y2 

d 

y2 

y3 

dt 

y3 

2sin(y4)-4y3-5y2-2y1 

1 

As  with  a single  equation,  we  approximate  the  dy  / dt  terms  using  1st  order  finite 
differences: 


yc"+1]  = yc'°  +h-  f (y*) . 

Accuracy,  Stability,  and  Step  Size 

Accuracy  is  a measure  of  how  closely  our  numerical  solution  y[n]  matches  our  “true” 
solution  of  y(t[n]  j . We  will  refer  to  two  types  of  accuracy.  The  first,  local  accuracy,  is  the 
accuracy  of  the  numerical  method  over  a single  time  step.  When  these  numerical  methods 
are  used,  we  will  take  multiple  time  steps.  The  errors  from  each  time  step  will  tend  to 
accumulate.  The  accuracy  of  a numerical  method  over  many  time  steps  will  be  referred  to 
as  the  global  accuracy.  We  will  sometimes  denote  this  global  error  as  sCn)  where: 

j/wSJ;(tCn))  + 8W  =>  s w=yw-y(tw) 
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where  y(tC/,))  is  the  true  solution  at  time  t[n] . 

Stability  is  a somewhat  ambiguous  term  and  may  appear  with  such  qualifiers  as  inherent, 
partial,  relative,  weak,  strong,  absolute,  etc.  In  general,  a solution  is  said  to  be  unstable  if 
errors  introduced  at  some  stage  of  the  calculations  are  propagated  and  magnified 
without  bound  throughout  the  subsequent  calculations.  An  unstable  result  will  not  only 
be  inaccurate  but  also  will  be  patently  absurd.  We  can  show  that  a criteria  for  stability  is: 

|sc"+1]| 


The  issues  of  accuracy  and  stability  are  important  when  trying  to  choose  one  or  more 
values  for  a step  size  for  integration.  We  will  generally  need  a step  size  much  smaller  than 
the  range  of  time  in  which  we  are  interested,  so  many  time  steps  must  be  performed  to  get 
our  numerical  answer.  The  step  size  chosen  should  be  small  enough  to  achieve  our 
desired  accuracy,  but  should  be  as  large  as  possible  to  avoid  an  excessive  number  of 
derivative  evaluations. 

In  general,  accuracy  relates  to  the  characteristics  of  the  method  when  making  the  step  size 
smaller  and  smaller,  whereas  stability  relates  to  the  characteristics  of  the  method  when 
making  the  step  size  larger  and  larger. 

Forward  Euler  Method 

The  simplest  method  is  to  evaluate  the  derivative  using  the  values  at  the  beginning  of  the 
time  step.  The  recursion  formula  is  then: 

yn+i) = yn) + h ■ f(/n)  ,t™)=yM+h-fM 

where  /(n)  refers  to  the  derivative  function  evaluated  using  the  n-th  time  step  values.  This 
method  is  often  referred  to  as  the  forward  Euler  method.  The  method  is  simple  and 
explicit  (i.e.,  requires  no  iteration  to  get  y[n+1] ).  Its  disadvantages  are  that  it  is  not  very 
accurate  (0(h)  global  accuracy)  & can  have  stability  problems  with  time  steps  that  are  too 
large. 

Let’s  look  at  numerically  solving  the  ODE: 

— = 1 - 2y  with  j/(0)  = 0 . 

Note  that  there  is  an  analytical  solution: 
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yW=i(i-e-“). 

The  forward  Euler  recursion  formula  will  be: 

yn+i]  = yn]  + ^i-2ynl)=/2+(i-2/2)_yc"->. 

The  following  chart  shows  the  numerically  approximated  solution  at  various  step  sizes  as 
compared  to  the  analytical  solution.  Notice  that  at  a reasonably  small  step  size  of  h = 0.1 
the  numerical  solution  is  pretty  good  — the  "knee"  of  the  curve  is  pretty  close  & we  get  to 
the  correct  ultimate  values.  However,  the  larger  step  sizes  give  unreasonable  values. 
h = 0.5  gets  to  the  ultimate  value  in  only  a single  step  & h = 0.9  oscillates.  In  fact,  h = 1 
oscillates  between  1 and  0 while  h > 1 oscillates  in  a diverging  manner. 

Forward  Euler  Solution  to  y'  - 1 - 2y,  y(0)  = 0 


The  problems  with  h > 1 show  the  stability  problems  of  the  forward  Euler  method.  Errors 
introduced  at  early  time  steps,  instead  of  dying  out,  start  to  grow  and  dominate  the 
numerical  solution.  This  shows  that  the  forward  Euler  method  is  conditionally  stable  for 
small  enough  step  sizes.  For  this  problem,  the  stability  limit  is  h = 1.  However,  what  if  we 
just  wait  until  the  end  of  the  problem  to  change  to  a large  step  size.  The  following  chart 
shows  what  happens  if  we  use  a step  size  of  h = 0.1  until  t = 2 & then  use  a step  size  of 
h = 1.1 . Notice  that  even  though  we  are  almost  at  the  ultimate  value,  the  large  step  size  is 
still  unstable. 
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Forward  Euler  Solution  to  y'  = 1 - 2y,  y(0)  = 0 
h = 0.1  for  f < 2 
h = 1.1  for  f >2 


t 


Backward  Euler  Method 

For  the  forward  Euler  method,  we  evaluated  the  derivative  function  at  the  beginning  of  the 
time  step,  /(y[n],tCn]).  But  we  can  evaluate  / at  any  appropriate  position.  Ifweusethe 
value  at  the  end  of  the  time  step  we  get: 

yn+1)  = yw  + h ■ f(y(n+1\t[n+1) ) = /"]  + h ■ f[n+1) . 

This  method  is  often  referred  to  as  the  backward  Euler  method.  The  method  is  still 
simple,  but  it  is  implicit  (i.e.,  it  may  require  iteration  to  get  y[n+1] ).  It  only  has  0(h)  global 
accuracy  just  like  the  forward  Euler  method,  but  it  has  much  better  stability  characteristics. 
For  our  example,  we  now  have: 

yM=yW  + h(  l-2/”+15)  =>  y^)=Zi!!l/7. 

v ' 1+2 h 

Note  that  in  this  particular  case  we  can  solve  directly  for  y[n+1] , but  this  is  not  always  the 
case.  The  following  chart  shows  the  numerical  solution  at  various  step  sizes  as  compared  to 
the  analytical  solution.  Notice  that  at  a reasonably  small  step  size  of  h = 0.1  the  numerical 
solution  is  pretty  good  just  like  for  the  backward  Euler  method  — the  "knee"  of  the  curve  is 
pretty  close  & we  get  to  the  correct  ultimate  values.  The  biggest  difference  is  that  even 
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though  larger  step  sizes  are  not  very  accurate,  the  results  do  not  oscillate  or  diverge.  It  can 
be  shown  that  the  backward  Euler  method  is  unconditionally  stable  for  any  step  size. 

Backward  Euler  Solution  to  y'  = 1 - 2 y,  y(0)  = 0 


Comparison  of  the  Stability  Characteristics  of  the  Euler  Methods 

A common  method  used  to  examine  the  stability  characteristics  of  a numerical  method  is  to 
determine  what  happens  when  the  method  is  applied  to  the  test  equation: 

^-  = -A y with  y(0)  = l 

where  X is  real  & positive.  This  test  equation  has  the  analytical  solution 
y(t)  = l-e~u. 

If  we  expression  the  solution  as  the  sum  of  the  exact  solution  and  the  error  s , then  the 
error  must  also  satisfy: 


dt 
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We  can  now  see  what  happens  to  the  error  along  a time  step.  For  the  forward  Euler 
method: 


8(n+D=eW+/7l- 


,(«+!)  _ eM 


h(-X  sCn)) 
(1  -Xh) 


So,  for  stability: 


,(»+!) 


,(«) 


= (l-Xh)\<l 


0 <h<- 
X 


Only  values  of  h<X/2  will  give  stable  results.  This  is  why  the  forward  Euler  method  is 

conditionally  stable. 

Remember  our  example  problem?  We  found  that  the  numerical  solution  oscillated  up  to 
h = 1 and  then  was  unstable  for  h > 1 . From  this  relationship  above,  the  example  problem 
must  have  had  A,  = 2 (remember  the  term  e 2t  which  corresponds  to  the  term  e~Xt  in  the 
test  problem?).  Also  note: 

g["+i)  =-ki-X/2)|st',)  when  — </?<  — . 

X X 

The  error  will  not  necessarily  grow,  but  it  will  change  sign,  leading  to  an  oscillatory 
behavior  that  is  not  present  in  the  exact  solution.  For  some  problems  these  oscillations 
may  be  noticeable  and  unacceptable. 

Now  lets  do  a stability  analysis  for  the  backward  Euler  method: 

8Cn+13=8t")+h(-XsC"+1)) 


(1  + Xh) 

So,  for  stability: 


g(n+l) 

1 

8Cn) 

(1  + Xh) 

h>  0 


Now,  all  values  of  h will  give  stability.  This  is  why  the  backward  Euler  method  is 
unconditionally  stable.  This  was  shown  in  our  example  problem. 
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Taylor  Series  Expansions  & Order  of  Accuracy 

We  can  use  Taylor  series  expansions  to  create  new  integration  schemes  and  determine 
their  order  of  accuracy.  Remember  that  the  Taylor  series  expansion  of  a function  around 

y(t) is: 

y(t  + h)  = y(t)  + hy\t)  + ^-y"(t)  + ^-y"\t)  + --- 

Z o 

We  could  also  have  written  the  expansion  around  y(t  + h ) : 

y{t)  = y(t  + h)  - hy'(t  + h)  + y'\t  + h) + h)  + - • . 

z 6 

Not  all  of  the  derivatives  of  y(t)  are  explicitly  known.  We  may  need  to  compute  the 
derivatives  needed  in  the  Taylor  series  expansion  from  the  defining  equation  itself.  Since: 

/(0=f=/(M) 


and: 


^ 5/,  a/, 

df  = —dy  +—dt 
oy  dt 


v = %&  + % = % f + W 

dt  dy  dt  dt  dy  dt 


then: 


y"(&- 

?"(*) 


^L=^LJlf+tL 

dt2  dt  dy  dt 


d^y_ 

dt 3 


rd2f  dy  d2f  A 


dy  dt 

eIff+ 


dydt 

d2f 


dy  dydt 


j dyydydt  dt  y 

fMf+tL)+ 

dy\dy  dt ) 


+ 


d2f  dy  + d2f 


dydt  dt  dt2 


dydt  dt 


\ 


and  so  on. 


Let  us  take  the  expansion  around  j/(t) . We  can  truncate  the  expression  after  the  linear 
term  and  get  the  forward  Euler  method: 

y(t  + h)«y(t)  + hy'(t). 
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Since  we  have  neglected  the  2nd  order  and  higher  terms,  it  is  said  that  the  approximate 
function  is  of  order  accuracy  0[h2) . This  represents  the  local  accuracy  (i.e.,  accuracy  along 
a single  step).  Generally,  the  error  starts  to  grow  when  more  than  one  step  is  made  — the 

global  accuracy  will  generally  be  one  order  less  than  the  local  accuracy.  So,  for  the 

forward  Euler  method,  the  global  accuracy  will  be  0[h ) . 

If  we  take  the  expansion  around  j/(t  + h)  and  truncate  the  expression  after  the  linear  term, 
we  get  the  backward  Euler  method.  Since  we  have  neglected  the  2nd  order  and  higher 
terms,  it  has  local  accuracy  of  order  0[h 2)  just  like  the  forward  Euler  method.  Also  like  the 
forward  Euler  method,  the  global  accuracy  will  be  0[h ) . 

Just  what  does  the  order  accuracy  mean  to  us?  If  a method  has  order  accuracy  0[hm) , then 
reducing  the  step  size  by  the  fraction  F will  reduce  the  error  in  the  numerical  solution 
from  s to  Fm s . For  example,  for  order  accuracy  0(h),  cutting  the  step  size  in  Vz  will  cut 
the  error  also  in  Vz.  However,  for  order  accuracy  0(h2),  cutting  the  step  size  in  V2  will  cut 
the  error  to  %.  This  shows  a real  advantage  of  methods  that  have  higher  orders  of 
accuracy. 

How  can  we  get  higher  orders  of  accuracy?  We  can  include  the  effects  of  the  higher  order 
derivatives  directly.  Or,  we  can  use  additional  function  evaluations  to  simulate  the  effects 
of  the  higher  order  derivatives. 

Runge-Kutta  Methods 

It  is  possible  to  develop  procedures  that  involve  on  1st  order  derivative  evaluations  but 
which  produce  results  equivalent  in  accuracy  to  the  higher-order  Taylor  expansion 
formulas.  These  algorithms  are  the  Runge-Kutta  methods.  Approximations  of  the  2nd,  3rd, 
and  4th  orders  are  equivalent  to  keeping  the  h2 , h3,and  h4  terms,  respectively.  These 
require  function  evaluations  at  2,  3,  and  4 points,  respectively,  along  the  interval 
t(n)  < t < ttn+1) . (Method  of  order  v > 4 require  more  then  v function  evaluations.) 

The  Runge-Kutta  formula  of  order  v can  be  written  as: 

= ym +£wlk,  ■ 

i= 1 


Where  the  w;.  are  weighting  factors  whose  sum  is  1 (i.e.,  w1  + w2  -1 — + wv  = 1 ) and: 


K=hf 


i- 1 


t[n]+cih,yw  + Yjai,Jki 

i= 1 J 
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and  where  the  cj  and  ai  . are  constants.  The  values  for  the  w, , cl , and  a,  . constants  are 
chosen  to  simulate  keeping  higher  order  derivative  terms  in  the  Taylor  series  expansion  & 
hence  give  high  orders  of  accuracy. 

The  Runge-Kutta  methods  belong  to  a class  of  predictor-corrector  methods.  The  first  step 
is  to  use  an  explicit  formula  to  predict  the  value  at  the  end  of  the  time  step.  Then,  implicit 
formulas  are  used  to  correct  this  initial  estimate.  However,  instead  of  iterating  to 
convergence  using  the  implicit  formula,  the  initial  estimate  is  used  for  only  a couple 
correction  steps. 

A 2nd  order  Runge-Kutta  method  is  known  as  Heun’s  method  or  the  improved  Euler 
method.  It  uses  a forward  Euler  prediction  step  with  a backward  Euler  correction  step.  The 
formula  is: 


/">=/'"+ \k+\k 

K=h-f{&\yw) 

k2=h-f(t^+h,yw+k1). 


The  stability  equation  for  this  method  is: 


£c-i) 


8 


(») 


1-—Xh 

2 


1 + -Xh 
2 


<1 


h>0 


but  will  lead  to  oscillations  for  h>  2/X. 

Perhaps  the  most  widely  used  Runge-Kutta  method  is  the  4th  order  method  with  constants 
developed  by  Gill.  The  method  is: 

(»  + D = yW  + K +-/f  +-  k. 3+~k, 

6 1 3 2 3 3 6 4 

K=h-f(tw,yM) 

k,=h-f[tw  + |A,y">  + i*r1\ 

k3  =h-fitw  +— h,yw  +akt  +bk2  . 

V 2 J 

k4  =h-  f(tw  +h,yM  +ck2  +dk3^j . 
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where  the  constants  are: 

V2-1  , 2-V2  V2  J 2 + V2 

a = , b = , c = , a = . 

2 2 2 2 

It  is  reported  that  the  Runge-Kutta-Gill  method  is  stable  for  h<  2. 8 /A,. 

Let’s  again  look  at  numerically  solving  our  sample  ODE: 

^-  = 1-2 y with  y(0)  = 0 . 

The  following  chart  shows  the  numerically  approximated  solution  using  the  2nd  order 
improved  Euler  method.  We  are  using  the  same  step  sizes  as  before.  Remember  that  for 
the  larger  step  sizes  in  the  forward  Euler  method,  the  numerical  solution  was  unstable. 
Here,  the  larger  step  sizes  are  not  very  accurate,  but  they  are  not  unstable.  Plus,  the 
smaller  step  size  of  h = 0.1  gives  a very  good  numerical  solution  — this  solution  is  visually 
better  than  either  of  the  forward  or  backward  Euler  solutions. 

Finally,  the  following  chart  shows  the  numerically  approximated  solution  using  the  4th 
order  Runge-Kutta-Gill  method.  Notice  that  the  numerical  solution  is  excellent,  even  for  the 
fairly  large  step  size  of  h = 0.5 . 

Runge-Kutta-Gill  Solution  to  y'  = 1 -2y,  y(0)  = 0 
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Stiff  Equations 

An  important  characteristic  of  a system  of  ordinary  differential  equations  is  whether  they 
are  stiff  or  not.  A system  is  stiff  if  the  characteristic  time  constants  for  the  problem  vary  by 
many  orders  of  magnitude.  Mathematically,  if 


dy 

dt 


f(y)*Ay 


by  linearization,  then  there  are  a set  of  eigenvalues.  A,(. , which  form  the  solution  to: 
det^A-X,.!  Wo 


The  solution  to  the  linearized  ODE  will  be  of  the  form: 


y< = Yjcie 


where  X.=  — . 

T- 


A stiffness  ratio  S can  be  defined  as: 

s _ max(Re(A,;.)) 
min(Re(?i;.)) 

Typically,  S > 1000  is  considered  a stiff  system  & S > 106  is  a very  stiff  system. 

What  does  this  mean  to  us?  If  a solution  is  made  up  of  a sum  of  exponential  decay  terms  all 
with  different  time  constants,  then  we  must  integrate  over  a time  length  that  will 
encompass  the  term  that  takes  the  longest  to  decay.  This  will  be  controlled  by  the  largest 
time  constant  (i.e.,  the  smallest  eigenvalue).  However,  for  stability  considerations,  the 
largest  step  size  that  can  be  used  is  controlled  by  the  smallest  time  constant  (i.e.,  the  largest 
eigenvalue).  When  these  are  different  by  several  orders  of  magnitude,  we  may  have  a 
problem  in  trying  to  perform  the  numerical  solution  in  a reasonable  number  of  time  steps. 

Lets  look  at  a stiff  example.  The  following  ODE  system: 


d 

y i 

-500.5 

499.5  “ 

y i 

where 

>.(»)' 

~2~ 

dt 

_y2J 

499.5 

-500.5 

_y2J 

1 

has  the  solution 
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"l.5e_t  + 0.5e~looot" 

l.Se1  -0.5e~looot_ 

Here,  the  eigenvalues  are  X = [l,1000]r  and  the  equivalent  time  constants  are 
t = [l,0.00l]  . The  stiffness  ratio  is  S = 1000,  signifying  a stiff  problem.  The  exponential 
term  associated  with  X = 1000  will  decay  99%  of  the  way  to  its  ultimate  value  (of  zero)  by 
t = 0.0046 ; however,  the  exponential  term  associated  with  X = 1 will  decay  99%  of  the  way 
to  its  ultimate  value  (of  zero)  by  t = 4.6 . If  we  use  the  forward  Euler  method,  we  must  keep 
the  step-size  as  h < 2 / A,max  to  keep  the  numerical  solution  stable.  For  this  problem,  that 
means  h < 0.002 , meaning  that  we  would  need  at  least  2,300  time  steps  to  reach  our  final 
time  of  t = 4.6 . 

The  following  figure  shows  the  numerical  solution  at  early  times.  In  this  chart,  the  step  size 
used  is  h = 0.0001,  much  smaller  than  the  stability  limit  of  h<  0.002.  However,  using  this 
extremely  small  step  size  will  require  46,000  steps  to  reach  the  ultimate  limit  of  t-  4.6 . 
What  happens  if  we  raise  the  step  size  after  t > 0.01 , where  the  effects  of  the  larger 
eigenvalue  has  already  died  out? 


Forward  Euler  Solution  to  Sample  Stiff  ODE  System 

/r=0.0001 


t 


The  next  figure  shows  the  numerical  solution  after  t > 0.01  when  the  step  size  is  increased. 
We  will  go  just  above  the  stability  limit  of  h<  0.002  to  h = 0.00201.  Notice  that  just  this 
little  bit  above  the  stability  limit  starts  to  show  instability  — it  takes  a while,  but  the 
oscillations  do  occur  after  t > 1.0 . 
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Forward  Euler  Solution  to  Sample  Stiff  ODE  System 
h =0.0001  for  t <0.01 
/r=0.00201  for  f>0.01 

2.5 

2.0 

1.5 

y(t) 

1.0 
0.5 
0.0 

0.0  0.2  0.4  0.6  0.8  1.0  1.2  1.4  1.6  1.8  2.0 

t 

So  how  do  we  numerically  solve  stiff  problems?  We  need  to  use  methods  that  are  of  high 
order  but  also  unconditionally  stable.  This  means  that,  in  general,  some  type  of  implicit 
method  will  have  to  be  used.  Methods  have  been  published  by  C.  W.  Gear  and  are  available 
as  the  GEARB  package  of  routines.  These  routines  have  proven  to  be  very  popular  for 
solving  stiff  problems. 

Some  of  the  newer  mthods  that  are  successful  in  solving  stiff  systems  are  based  upon 
lineaizing  the  set  of  derivatives  and  effectively  doing  a single  Newton  iteration  at  the  end  of 
the  step.  These  are  referred  to  as  semi-implicit  methods.  The  following  system  of 
equations  is  linearized: 

y(n+l)  =yW  + /jf  (y["+1^  =^>  yt"+13  ~yC»)  + /2  f^yC'0^+^l  (y(n+0  _ yW ) 

3y  y(») 

and  then  can  be  rearranged: 

I-/2—  [~ y [n+1)  _ y (»)  ] = h f ( y C'°  ) 

8y^\L  J v ' 

and  solved  as: 
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— I"1 

y("+1)=yw+/i  I-/3—  f(yw) 

v ' 

Here,  the  derivative  term  is  the  Jacobian  matrix,  J , with  elements  ]tj  = dfj  / <3y . . The  final 
solution  for  yc"+1)  requires  some  type  of  matrix  inversion  (usually  an  effective  matrix 
inversion  via  an  LU  decomposition). 

Rosenbrock  methods  are  semi-implicit  Runge-Kutta  methods.  One  particular  set  of 
methods  have  been  proposed  by  Kaps  & Rentrop.  Their  recursion  equations  have  the 
following  form: 

4T-F>  g,  = f(yw 

y h J v 

^ T-Tw  + 

— I - JCn)  g3  = f (yC"}  + a3lSl  + «32g2  ) + ^[c3lgl  + C32g2  ] 

— I - JCn)  g4  = f (yw  + a31gx  + a32g2 ) + -[c41g4  + c42g2  + c43g3 ] 

yCn+1]  = y(n)  + ^igi  + b2%i  + kSs  + KZa 

where  J(n)  is  the  Jacobian  matrix  evaluated  at  the  beginning  of  the  step. 

The  following  Table  1 shows  the  parameters  for  this  method.  Note  that  there  are  two  sets 
of  parameters,  one  set  for  a 4th  order  method  and  another  set  for  a 3rd  order  method. 
These  two  sets  of  methods  are  very  useful  since  they  can  be  used  to  estimate  errors  on  a 
particular  step  as: 

77  _ lirtn+l)  _ vtn+1)| 

c \y  4th  J3rd  \ 

~ |(^l,4tfc  _ ^l,3rd)gl  +(^2,4 th  ~ ^2,3rd)g2  (^3,4tft  _ ^3,3rd)g3  +(^4,4t/i  ~^4,3rd)g4  — O^h  ) 

and,  hence,  lead  to  an  algorithm  for  increasing  or  decreasing  the  step  size  to  achieve  a 
particular  tolerance. 
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Table  1. 

Kaps-Rentrop  Semi-Implicit  Runge-Kutta  I 

Parameters 

Parameter 

Value 

Parameter 

4th  Order 
Value 

3rd  Order 
Value 

a21 

2 

\ 

19/9 

97/54 

a31 

1.92 

b2 

1/2 

11/36 

a32 

0.24 

\ 

25/108 

25/108 

C21 

-8 

K 

125/108 

0 

C31 

14.88 

C32 

2.4 

C41 

-0.896 

C42 

-0.432 

C43 

-0.4 

Y 

0.5 
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Summary  of  Formulas 

Forward  Euler  Method: 

/n+1)  = yw  + h ■ f(yln\tw ) = y[n]  + h ■ f[n) 

Backward  Euler  Method: 

yn+i) = yn) + h ■ / (yn+1]  ,tCn+l) ) = yn) + h ■ /Cn+1) . 

Improved  Euler  Method: 

y",=/",+y+y 

where:  k1=h- f(t[n]  ,y[n)^J 

k2=h-f(tw+h,yw  + k1) 

Runge-Kutta-Gill  Method: 

/“,=/">+y+y+y+y 

0 5 5 0 

where:  k1=h- f(t[n]  ,y[n)^J 

k2=h-f(tM  + h,yM+ht\ 

k3=h- f(tw +—h,yw +ak1+bk2  . 

V 2 y 

k4  =h-  f(t[n]  +h,y[n)  +ck2  +dk3^j . 

V2-1  , 2-V2  V2  J 2 + V2 

a = . b = , c = , d = . 
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